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Modelling the flow in a thin accretion disc like a dynamical system, we analyse the nature of 
the critical points of the steady solutions of the flow. For the simple inviscid disc there are two 
critical points, with the outer one being a saddle point and the inner one a centre type point. For 
the weakly viscous disc, there are four possible critical points, of which the outermost is a saddle 
point, while the next inner one is a spiral. Coupling the nature of the critical points with the outer 
' boundary condition of the flow, gives us a complete understanding of all the important physical 

features of the flow solutions in the subsonic regions of the disc. In the inviscid disc, the physical 
realisability of the transonic solution passing through the saddle point is addressed by considering 
a temporal evolution of the flow, which is a very likely non-perturbative mechanism for selecting 
the transonic inflow solution from among a host of other possible stable solutions. For the weakly 
viscous disc, while a linearised time-dependent perturbation imposed on the steady mass inflow rate 
causes instability, the same perturbative analysis reveals that for the inviscid disc, there is a very 
close correspondence between the equation for the propagation of the perturbation and the metric 
oo : of an acoustic black hole. Compatible with the transport of angular momentum to the outer regions 

of the disc, a viscosity-limited length scale is denned for the spatial extent of the inward rotational 
drift of matter. 
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In astrophysical fluid dynamics, studies in accretion processes occupy a position of prominence. Of such processes, 
critical (transonic) flows — flows which are regular through a critical point — are of great importance . A classic 
example of such a flow is the Bondi solution in steady spherically symmetric accretion 0, 0. What is striking 
about the Bondi solution is that while the question of its realisability is susceptible to great instabilities arising from 
■ infinitesimal deviations from an absolutely precise prescription of the boundary condition in the steady limit of the 
hydrodynamic flow, it is easy to lock on to the Bondi solution when the temporal evolution of the flow is followed. 
This dynamic and non-perturbative selection mechanism of the transonic solution agrees quite closely with Bondi's 
conjecture that it is the criterion of minimum energy that will make the transonic solution the favoured one 

The appeal of the spherically symmetric flow, however, is limited by its not accounting for the fact that in a realistic 
situation, the infalling matter would be in possession of angular momentum — and hence the process of infall should 
lead to the formation of what is known as an accretion disc. A simple and well-understood model of such an accreting 
system is the inviscid and thin accretion disc 0, 0, El 0> 13 j an d this we will initially consider in our study to provide 
ourselves with an effective guideline for a subsequent approaching of the more general viscous problem. However, the 
inviscid model has its own limitations because, while the presence of angular momentum leads to the formation of 
an accretion disc, a physical mechanism must be found for the outward transport of the angular momentum, which 
should then make possible the infall of the accreting matter into the potential well of the accretor. Viscosity has been 
known to be such a physical means to effect infall, although the exact presciption for viscosity in an accretion disc 
is still a matter of much debate What is well appreciated, however, is that the viscous prescription should be 
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compatible with an enhanced outward transport of angular momentum. The very well known a parametrisation of 
Shakura and Sunyaev [9|, is based on this principle, and we invoke this a model to logically extend our study from 
an inviscid disc to a "quasi-viscous" disc. With that objective we have prescribed a first-order viscous perturbation 
(dependent on a) about the zeroth-order inviscid solution, and we see that the merest presence of viscosity in the 
flow equations, gives rise to significant qualitative changes in the flow behaviour, insight about which could not have 
been derived from the inviscid disc model alone. 

In our present study, the accretion flow has been considered to be driven by the classical Newtonian potential. 
For accretion on to an ordinary star or even a white dwarf, this is arguably a satisfactory working premise. We 
set up the governing steady flow equations of both the inviscid and the quasi-viscous disc like a dynamical system. 
This approach gives us a clear understanding of the nature of the critical points of the flow, as well as any sense of 
direction that may be associated with each solution. Combining this understanding with a knowledge of the outer 
boundary condition for physically meaningful inflow solutions, makes it possible to construct the phase diagram of 
the flow solutions for both the cases we are studying. For the inviscid disc it has been found that there are exactly 
two critical points for the steady flow solutions, of which, for realistic boundary conditions, the outer one is a saddle 
point, while the inner one is a centre type point. The situation gets somewhat more complicated for the quasi- viscous 
disc. In this second case we see that there should be four critical points, of which we analyse the behaviour of those 
two which, in the limit of vanishing a would take us back to the phase portrait of the inviscid flow. The outer one 
of these two critical points remains saddle type as before, while the inner one behaves like a stable spiral. While 
dwelling on this matter, we should mention that studies on the nature of the critical points in an accretion disc have 
been reported before. While numerically studying the effect of a-viscosity on the flow structure around the inner 
edge of the disc near a non-rotating black hole, Matsumoto et al. discussed the behaviour of critical points for 
the phase portrait of the isothermal flow, and for small a they upheld the case for a transonic saddle point. Along 
similar lines Afshordi and Paczyhski [T^ adopted a semi-analytic approach to understanding the nature of the fixed 
points for constructing a physical solution. In making a quantitative analysis of the inner boundary condition, for a 
constant effective speed of sound, they also found the existence of a saddle point for small a. While mentioning these 
background issues, we must stress however, that our interest in the disc would be to study its stability aspects and 
its long-time evolutionary properties on large length scales, and to do so we would also need to have a clear notion of 
at least the qualitative features of the phase portrait of the steady flow solutions, and its critical points. That these 
important and qualitatively physical conclusions about the flow could be drawn without taking any recourse to the 
conventional practice of a numerical integration of the steady flow equations, amply demonstrates the simplicity, the 
elegance and the power of the dynamical systems approach that we have adopted to study the thin accretion disc. 

Following our understanding of the nature of the critical points in the phase portrait, the question that we then take 
up is about the preference of the accreting system for any particular velocity profile in the stationary phase portrait, 
and a selection criterion thereof. In this regard the transonic solution has always been the favoured candidate. 
However, it is common knowledge from the study of dynamical syste ms, that a flow solution (the transonic flow in 
our case) passing through a saddle point cannot be realised physically [13j . To address this issue satisfactorily it must 
be appreciated that the real physical flow is not static in nature, but has an explicit time-dependence. With respect 
to this point it is tempting to subject the steady flow solutions to small perturbations in real time, and then study 
their behaviour. This has been done elsewhere |l4| for the inviscid disc, and it has been shown that the steady inflow 
solutions of abiding interest are all stable under the influence of a linearised time-dependent perturbation on the mass 
inflow rate. The application of a similar treatment on the quasi-viscous disc, however, very much differently affects 
the stability of its steady inflow solutions, which, as we have shown here, become destabilised by small perturbations 
— treated both as standing and travelling waves. 

Since one way or the other, no direct conclusion could be drawn about the selection of a particular solution through 
a perturbative technique, for the inviscid disc at least, we try to have an understanding of a true selection mechanism 
and the attendant choice of a particular solution, by studying the evolution of the accreting system through real time. 
A model analog shows that it is indeed possible in principle for the temporal evolution to allow for the selection of an 
inflow solution that passes through the saddle point, and under restricted conditions we demonstrate that the selection 
criterion conforms to Bondi's minimum energy argument, which is invoked to favour the transonic solution in the 
spherically symmetric case. Interestingly enough in this context, we have also shown that although the perturbative 
study has offered no direct clue about the selection of a solution, the equation governing the propagation of an acoustic 
disturbance in the flow bears a close resemblance with the effective metric of an acoustic black hole. We make use of 
this similarity to argue that the flow would cross the acoustic horizon transonically. 

For the quasi-viscous disc, dissipation upsets all idealised conditions, and allows for the setting of no precise 
criterion for selecting any particular solution, something that we might otherwise establish under conserved conditions. 
Nevertheless, we have at least been able to argue that with the cumulative transfer of angular momentum to the outer 
regions of the disc, a characteristic scale of length may be set, beyond which the accumulation of angular momentum 
resists further inward drift of matter. 
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II. THE STEADY EQUATIONS OF THE FLOW AND ITS CRITICAL POINTS 

The flow that we consider is an axisymmetric thin disc whose local height is given by H |9j. It is customary to 
describe the steady flow by the equation of continuity, the equation for radial momentum balance and the equation 
for angular momentum balance, all of which are respectively given in the following notation 

by 
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In the above, v is the radial drift velocity, R is the radial distance, Q is the angular velocity, 17k is the Keplerian 
angular velocity defined as = GM/R 3 , c s is the velocity of sound defined as c 2 = dP/dp and a is the effective 
viscosity parameter of Shakura and Sunyaev @, • 

It is a standard practice to make use of a general polytropic equation of state P = kp 1 where k and 7 are constants, 
with 7 being the polytropic exponent, whose admissible range (1 < 7 < 5/3) is restricted by the isothermal limit and 
the adiabatic limit, respectively. The condition of hydrostatic equilibrium along a direction perpendicular to the thin 
disc, allows us to use the approximation [|J 



where vk — RQk- In that case Eq.([Q) leads to 
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in which n = (7 — 1) . The constant of integration in Eq.JSJ) can be physically identified with the mass inflow rate. 
Combining Eqs.Q and JSJ gives, 



djy)_ 

dR 



2v^ 



5c s 2 -( 7 + l)(^-tt 2 ) R 2 



(7 + 1)' 



2c? 



and a first integral of Eq. Q can be written as 
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in which L is a constant of integration. 

The inviscid limit of Eq.Q may be obtained by setting a — 0. This model has found some use in accretion 
studies 0,0, and it gives us the condition flR 2 = L, where L, which can be physically identified as the specific 
angular momentum of the flow, now becomes a constant of the motion. Such a constraint gives us an explicit 
dependence of f2 on R, and allows us to fix the critical points of the flow in the v 2 — R space. At such points both 
the numerator and the denominator in the right hand side of Eq. © vanish simultaneously Q , and this will deliver 
the critical point conditions as 
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with the subscripted label c indicating critical point values. The first condition gives the flow velocity at the critical 
point in terms of the local speed of sound at the same point. The latter condition, when written down explicitly as 
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gives a dependence of R c on c sc . This is a quadratic in R c , indicating that there are two critical points. The critical 
points (v 2 ,R c ) could be fixed in the v 2 — R space, if c sc could be expressed in terms of the constants of the flow 
system. To do so, it should be necessary for us to look at the integral of Eq.(J2J), which reads as 

'l + JL-2M. + j^ = E (10) 

2 2R? R 7 - 1 v ' 

in which the constant of integration E, which is actually the Bernoulli constant, can be determined with the help of 
the boundary condition that for very great radial distances, the speed of sound approaches a constant "ambient" value 
c s (oo), while the radial drift velocity v becomes vanishingly small. On using the critical point conditions in Ea. (|10|l . 
it could be seen that c sc could in principle be represented solely in terms of the physical constants of the system. 
Futhcrmore, for this flow system, we require that c s should be a monotonically increasing function for decreasing R, 
except perhaps very close to the accretor 0, and therefore c s should be single- valued at all R. On physical grounds 
alone, this should be a perfectly valid expectation. This would imply that for every value of i? c , there should be a 
single physical value of c sc , given by the profile of c s , and this value of c sc , having been fixed in terms of the system 
constants, would consequently enable the critical points to be fixed as well in the phase space of the flow solutions. 
This conclusion may be alternatively and equivalently derived by fixing c sc with the help of Eq. (JSJl . 
Thus the two solutions of the quadratic in R c , as given by Eq.JSJ), are obtained as 



_ 7 + 1 GM 
c ~ 10 cl 
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and the two fixed points may be considered to be at radii R = R c \ and R = R C 2, with R C 2 > R c i- We note from 
Ea. (|ll(l that real values of R c would necessitate the condition that for a regular solution passing through a critical 
point, there would be an upper bound to the admissible value of the constant angular momentum L. 

So far we have used the inviscid solution Q = LR~ 2 , derived by requiring that a = 0. About this inviscid solution 
for small non-zero values of a, we now prescribe a "quasi-viscous" solution, whose form, to a first order in a, we may 
write as 

n=-^+afl (12) 
R z 

in which the unknown function f2 is to be determined by invoking the angular momentum balance condition. To that 
end we use Ea. (|12fl in Eq.0, and by dropping all orders in a higher than the first, we will get 

vvkR 2 ^ 
which will allow us to define an effective angular momentum as 

, 2 „ „ A 2ac 2 
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Using this new functional dependence of Q in the radial momentum balance condition, given by Eq.®, we can account 
for the manner in which outward angular momentum transport effects the inward drift of matter. The validity of 
this dependence may be put to the test by studying the asymptotic relation of the angular momentum drift. On very 
large length scales the drift velocity v goes asymptotically as R~ 5 / 2 . Bearing in mind that for inflows, v is negative, 
on scales of length where the fluid assumes ambient conditions, we will have 

R- 3 



L cS - L + 2aL (15) 

Here Rl is a suitable scale of length, which, to an order-of-magnitude is given by R\ ~ G 'M m[Cg (oo) p(oo)] -1 , with 
Til being the mass inflow rate. This asymptotic behaviour is what we should expect entirely, because the physical role 
of viscosity is to transport angular momentum to large length scales of the accretion disc. 

Going back now to the critical point conditions, given by Eq.JSJ), we see that while the form of the first condition 
remains unaltered, the second condition, under the quasi- viscous regime (i.e. first order in a), gives 
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in which vk c = V GM Rc 1 is the Kepler velocity at the fixed points. We can therefore, eventually write Ea. l)lfi|l as a 
fourth-degree equation in R Cl provided, as before for the inviscid case, c sc is fixed in terms of the system constants. 
To do so, unlike before, we can only turn to the continuity condition, given by Eq.JSJ, since with viscous dissipation 
now present, the radial momentum balance equation cannot be exactly integrated. With c sc thus fixed, we conclude 
that in this new quasi-viscous scenario, there should be four critical points. Of these, we would be interested in those 
two points which, in the limit of a — > 0, would satisfactorily reproduce all the features of the inviscid flow. 



III. THE ACCRETION DISC AS A DYNAMICAL SYSTEM 



An integration of the first-order differential equation, given by Eq.©, for both the inviscid and the quasi- viscous 
disc, should in principle give us all the possible flow solutions in the v 2 — R space, corresponding to every possible 
value of the integration constant (itself to be determined with the help of boundary conditions for physical flows). 
This, however, could only be done numerically because the speed of sound is a function of both the radial drift velocity 
and the radial distance. So far this has been the conventional method of studying this problem. We contend here that 
alternatively a complete understanding of the behaviour of the flow solutions could be had, if standard techniques of 
dynamical systems analysis could be introduced to study an inviscid and thin accretion disc. To do that — having 
first determined the number of critical points — it would then be important for us to identify the nature of each of 
the critical points of the flow system. To that end we turn to Eq.© and parametrise it to get 
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This parametrisation has been carried out in a mathematical parameter space r, and in such a space the two 
equations above represent an autonomous first-order dynamical system \\',$ , given in the standard mathematical form 
y = Y(x, y) and x = X(x, y). To analyse the nature of the fixed points in this parameter space, v 2 and R would have 
to be expanded and then linearised in the expanded terms in Ea.l|17|). about the fixed point coordinates (R c ,v 2 ). 

For the inviscid disc, in terms of the expanded quantities Sv 2 and SR, a set of linear equations will then be given 
by 
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For the pair of linear first-order differential equations above, we use solutions of the form e for the perturbed 
quantities and this would deliver the eigenvalues A — growth rates of Sv 2 and SR in r space — as 



X 2_ A 5-7 GM n2 
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c- 
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where £ = (5 — 37) (5 — 7) _1 and Lkc = V GMR C is the local Keplerian angular momentum at the fixed points. 

The term in the square brackets in Ea. H19|) will determine the sign of A 2 . Knowing that two adjacent fixed points 
cannot be of the same nature fl3T ]. i.e. they will differ in their respective signs of A 2 , it is now evident that if R c is 
the outer fixed point R C 2, then A 2 will be positive and the fixed point will be saddle type, while if R c is the inner 
fixed point i? c i, then A 2 will be negative and the fixed point will be centre type. The flows of physical interest here 
are governed by the boundary condition that at large distances the drift velocity becomes vanishingly small, while 
the speed of sound approaches a constant value. This knowledge of the boundary condition, in conjunction with the 
nature of the fixed points, makes it possible to draw the phase trajectories of the flow, which are as shown in Fig. ^ 
The arrows in the figure, indicating the direction that may be associated with each flow solution, could only be drawn 
by taking recourse to a dynamical systems analysis. 

In a manner of speaking, the fixed point C2 at R = i? C 2, which is a saddle point, may be dubbed the "sonic point" 
because one of the two trajectories passing through it is a transonic inflow solution, rising from subsonic values far 
away from C2 to attain supersonic values for R < i? c2 . Here transonicity is to be attained when the drift velocit 
equals the speed of the propagation of an acoustic disturbance, which, for this system is given by \/2(j+ l) _1 / 2 c s 
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FIG. 1: A schematic diagram of the accretion problem. The fixed point C2 is a saddle point. Ci is a centre type point. The 
separatrices pass through C2. 



The transonic trajectory is shown as the heavy solid curve, and is an accretion inflow solution. Once the flow is in the 
supersonic region, the drift velocity starts decreasing for R < R c i. This happens because the centrifugal force in the 
flow (connected to its angular momentum, which is not lost under inviscid conditions) builds up resistance against 
the attractive influence of gravity. 

We can now study how these idealised conditions are affected when viscosity is brought into play. This we can do 
by using Eq. (|14|l to substitute for £1 in Eq. (|17|l . In carrying out a first-order perturbative expansion about the critical 
point values, as we have done before for the inviscid case, we will get the relation 
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The form of the second parametrised relation from Eq. i|17|) remains unchanged even with the introduction of viscosity, 
and this, together with Ea. i|2U|) . will eventually deliver a quadratic equation for the eigenvalues that reads as 
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We can now solve the quadratic in A to write, 



A = aX 1 ± 



Xo + aXo 



(23) 



in which, the viscosity parameter a can be tuned to arbitrarily small (though not exactly zero) values. This should be 
consistent with our quasi-viscous approximation, and this will also make Xi the dominant term in the discriminant 
of Ea. H23|l . As before, we can make use of the fact that no two adjacent critical points can be of the same nature. 
Combining this with the standard boundary condition, we can once again say that the outermost critical point will 
be a saddle point. However, an interesting change appears in the next inner critical point. For a — 0, we have seen 
that this point will be centre type, as has been shown in Fig. ^ Now a look at Ea. l|23|l immediately reveals that for 
very small non-zero values of a, with X\ being real and X2 < 0, this point becomes a spiral. This behaviour has been 
diagrammatically represented in Fig. [21 an d one can intuite that this is physically what it should be. With viscosity 
present, there is dissipation in the system, and the solutions therefore spiral in, as opposed to the conserved inviscid 



FIG. 2: A schematic diagram of the of the flow solutions in a weakly viscous accretion disc. The outer fixed point is a saddle 
point, but the next inner fixed point behaves like a spiral in the presence of viscous dissipation. 



situation, where the solutions about the point Ci form closed paths. Another important point to note is that with 
dissipation being present in the system, there would have been no general means for the integration of the equation 
for momentum balance. In that event, this dynamical systems analysis affords an easy way of studying the system. 

The dynamical systems approach, however, gives rise to a serious difficulty. The demonstration that the outermost 
fixed point is a saddle point, concomitantly places great difficulties in the way of realising the transonic trajectory 
passing through this point. This fact should be immediately evident from the direction that the arrows (characteristic 
of a saddle point) indicate for the solutions passing through C2 in Fig. ^ Conventional wisdom about the nature of 
saddle points leads to the understanding here that in the steady limit of the physical flows, the transonic solution 
would not be realisable [l3j . This issue may be understood in two ways. The first is more physical and is related 
to the prescription of the boundary condition. A look at Fig. ^ makes it quite obvious that the transonic solution 
could only be achieved after an infinitely precise determination of the outer boundary condition needed to generate it. 
Any deviation, however infinitesimally small, would take the system far away from transonicity. In the astrophysical 
context, it may easily be imagined that it is well nigh impossible to have such precise fulfilment of the boundary 
condition for transonicity. 

The second condition against achieving transonicity is somewhat more mathematical in nature. To find out the 
eigenvalues of the set of linear equations given by Ea . (|18|) . we have used a solution of the form SR ~ e Ar (the same 
form being used for Sv 2 as well), which we can recast as |r| ~ A -1 In \6R/R\, in which R is a suitable length scale 
factor to render the argument of the logarithm dimensionless. In the Sv 2 — SR space, reaching the saddle point (0, 0), 
would then entail an infinitely great divergence of |t|. This shows that any attempt to numerically integrate Eq.© 
would be thwarted even if the practical impossibility of determining an absolutely precise outer boundary condition 
were to be achieved. In Fig. this fact is represented by the way the arrows behave along the critical solutions, 
establishing that they are not physical solutions per se, but are actually separatrices of various classes of solutions, 
and this special status of the separatrices, as distinct from all other solutions, has been emphasised by the heavy solid 
curve and the heavy dashed curve passing through C2 . These arguments evidently apply for the saddle point even in 
the presence of viscosity. 

However, it is presently a firmly established fact that transonic solutions do occur in astrophysical systems such as 
the one being studied here. The difficulty about the realisability of the transonic solutions could be suitably addressed 
by appreciating that the discussion presented so far has been confined within the steady (time-independent) framework 
of the hydrodynamic flow, and that the arrows representing the direction of the flows are in a mathematical parameter 
space. On the other hand, in the real astrophysical context, the accreting system has an explicit time-dependence, and 
is evolving through time. Taking this fact into account could therefore satisfactorily resolve all difficulties regarding 
attaining transonicity. 

In passing, we make a point of academic interest. The nature of the outer fixed point C2 in Fig.^ has a bearing 
on the possible range of values that the constant specific angular momentum L may be allowed. For solutions passing 
through the outer critical point, we can easily recognise from Eq.(|SJ that the flow would be sub-Keplerian I n 
that case, we see that the condition (£/Lrc2) 2 < 1 would hold good. In addition to this, from Ea.l(T9)l. the saddle type 
behaviour of the point C2 would also imply that there would be another upper bound on L, given by (L/Lk C 2) 2 < C 
For the admissible range of the polytropic index 7, the possible range for £ would be < C < 1/2- This would 
then imply that the latter bound on L would be more restrictive, as compared to the former. It is interesting that 
this essentially physical conclusion could be drawn from modelling the disc as a dynamical system in a mathematical 
parameter space. 
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IV. A LINEARISED PERTURB ATIVE ANALYSIS IN REAL TIME 



As a preliminary study, the issue of explicit time-dependence in the accretion disc could be taken up from the 
viewpoint of a perturbative analysis in real time. In the thin disc approximation, with the help of a polytropic 
equation of state, the time-dependent generalisation of the radial momentum balance condition can be rendered as 



dv dv 2 dp GM 



L " s =0 



while the corresponding generalisation of the continuity equation can be set down as Rdt(pH) 
the latter, in terms of a new variable defined as f — p^ +1 '' 2 vR 5 ' 2 , can be further recast as 
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The angular momentum balance condition may similarly be represented as 
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(pvRH)' = 0, and 
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(26) 



From Eq. (|25|l it is obvious that the steady solution of /, given as /o, would be a constant, which might, within a 
numerical factor, be physically identified with the steady mass inflow rate. For orbits in a fixed gravitational potential, 
the smallness of the change in angular velocity allows us to ignore its explicit time variation [9j, [l6| ■ With the use of 
our quasi-viscous prescription, R 2 VL = L + aR 2 il' , in Eq.JSBJ), we then get 
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where L' is an integration constant. The effective angular momentum is now defined as L e g = L + aR 2 fl' and this is 
what we use in Ea. (|24|) after linearising in a. It is easy to see that in the steady limit Q' should assume the form of 
f2 as given by Eq.(JT3J 

The steady state solutions are given as vq, po and /o, and about these we impose small time-dependent perturbations 
v' , p' and /'. In terms of /', we can then derive an equation for the perturbation that is given by 



dt 2 



'dR 



df_ 

dt 



Vq dR 

4aL 2 



^0 (vq 



15 Cs0 > dR 



da 
dR 



37- 1 
7+1 



vo 



df_ 

dR 







(28) 



in which a = c^ /{vqvk) and (5 2 = 2/(7+1). 

We treat the perturbation as a standing wave and confine our analysis to the outer subsonic region of the flow. 
From Fig. [21 one can see that in the vicinity of the focus the solutions would be double- valued. This is physically not 
feasible, and therefore for all the solutions which spiral in towards the focus, their inner branch (for R < R c \) and 
their outer branch (for R > i? cl ) would have to be fitted via a vertical discontinuity in the neighbourhood of the focus. 
One may conceive of this discontinuity as a standing shock. For R > R c i, we will have a family of continuous solutions 
in the subsonic region, extending from the shock to very great distances from the accretor. Somewhere in this range, 
at two chosen points, we can constrain the perturbation spatially by requiring the standing wave to die out at those 
points. The outer point could suitably be chosen to be at the outer boundary of the flow itself, where by virtue of 
the boundary condition on the steady flow, the perturbation would naturally decay out. We choose the inner point 
to be infinitesimally close to the shock front, which can behave like a fully absorbing wall for all disturbances and 
make the perturbation vanish here too. Between these two points, the time-dependent behaviour of the amplitude of 
the standing wave will be an indicator of the stability of the steady solutions. For the perturbation we use a solution 
of the kind /' = p(R) exp(— iuot). Subsequently we multiply the result that we get from Eq.J2BJ by fog, and then 
carry out an integration by parts. This, followed by requiring that all the integrated "surface" terms vanish at the 
two boundaries of the standing wave, will give us a quadratic in lo, and it will be easy to see that —iuj will have a 
viscosity-dependent real part, even for arbitrarily small non-zero values of a. This real part is given by 



SJ(-iw) = 2a[ vopMR 



v p 



dR 



(29) 
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If 5ft(— iui) > 0, then the standing wave will have an amplitude growing in time. The term in the square brackets 
above will be crucial in this regard, and from its long-distance asymptotic behaviour at least, it is seen to be positive. 
Therefore on large length scales this indicates instability. 

To bolster our argument for this instability, we present the salient features of an analysis in which the perturbation 
has been treated as a radially travelling wave. To that end we prescribe the spatial part of the perturbation as 
p{R) — e s , where s is given by the power series s — J2^=-i k n (R)uj~ n . We confine our analysis to the short 
wavelength regime, in which the wavelengths are considered to be small compared to a relevant length scale in the 
system. The radius of the accretor, i?*, could be chosen as such a limiting length scale. 

To make any further progress we need to examine the structure of the function s. For that we go back to Eq.JSSJ) 
and use the solution f'(R,t) — exp(s — iujt) in this equation. In the resulting power series expansion, we sum up all 
the coefficients of the same power in u>, and set each individual sum to zero. In that way all the coefficients involving 
lo 2 will deliver the result 



v ± (3c s0 



dR 



while the coefficients of lo will lead to 



fco = -- In (f3v c s0 ) 



± aL 2 (3 7 - 1) 



(3c s0 {vp ± Peso) 
v v K R 3 Oq - /3 2c so) 



dR 



(30) 



(31) 



The two expressions above give the leading terms in the power series of s. For self-consistency it will be necessary 
to show that successive terms follow the condition oj~ n \k n {R)\ Lo~( n+1 ^\k n+ i(R)\. In the inviscid limit, this 
requirement can be shown to be very much true asymptotically for the first three terms, because we get fc_i ~ R, 
ko ~ In R and k% ~ R -1 . With the inclusion of viscosity as a physical effect, we see from Eas. (|30|) and l|31|l respectively, 
that while remains unaffected, ko acquires a viscosity-dependent term that goes asymptotically as aR. This in 
itself is an indication of the extent to which viscosity might alter the inviscid conditions. However, since we have 
chosen a to be very much less than unity, and since we have also chosen only short wavelength perturbations, implying 
that to 3> (vq ± /3c s o)/i?*, our self-consistency requirement still holds. Therefore, in so far as we are looking for a 
qualitative understanding of the effect of viscosity, it should be sufficient for us to consider the two leading terms 
only in the power series expansion of s, and with the help of these two, one may then set down an expression for the 
perturbation as 



f(R,t) 



A, 



e~ tuJt exp 



V ± PCsQ 



± aL 2 (37 - 1) 



/3c s0 (vp ± f3c s0 ) 
v v K R 3 (v 2 , - f3 2 c 2 ) 



dR 



(32) 



which should be seen as a linear superposition of two waves with arbitrary constants A + and A—. Both these two 
waves move with a velocity c s o relative to the fluid, one against the bulk flow and the other along with it, while the 
bulk flow itself has a velocity vq . It should be immediately evident to us that all questions pertaining to the growth 
or decay in the amplitude of the perturbation will be crucially decided by the real terms delivered to us from ko. 
The viscosity-dependent term is especially influential in this regard. For the choice of the lower sign in the real part 
of /' in Ea. (|32|l . i.e. for the outgoing mode of the travelling wave solution, we see that the presence of viscosity 
causes the amplitude of the perturbation to diverge exponentially on large length scales (where c s o ~ c s (oo) and 
since — Vq is positive for inflows. The inwardly travelling mode also displays similar behaviour, albeit 
to a quantitatively lesser degree. It is an easy exercise to see that stability in the system would be restored for the 
limit of a = 0, and this particular issue has been discussed elsewhere The exponential growth behaviour of the 
perturbation, therefore, is exclusively linked to the presence of viscosity. Regarding this point, we may also mention 
that Kato et al. [17j pointed out the existence of a non-propagating growing perturbation localised at the critical 
point, but interestingly enough they also found that this instability disappears in inviscid transonic flows. Chen and 
Taam ^| also made a numerical study of instability in the disc. For short wavelength radial perturbations, they 
found that the inertial-acoustic modes are locally unstable throughout the disc, with the outward travelling modes 
growing faster than the inward travelling modes in most regions of the disc, something that is very much in keeping 
with what Eq. l|32[l indicates here. 

Since we are more concerned with the question of the preference of the thin disc system for any particular solution, 
it now suffices for us to appreciate that in this regard, a perturbative time-dependent analysis offers no clue. For the 
inviscid disc all physically meaningful inflow solutions are stable, while with the introduction of viscous dissipation 
these solutions are destabilised on large length scales — a region through which all inflows will have to pass anyway. 

In concluding this section it would be instructive to furnish a parallel instance of the destabilising influence of viscous 
dissipation in a system undergoing rotation : that of the effect of viscous dissipation in a Maclaurin spheroid |19| . 
In studying ellipsoidal figures of equilibrium, Chandrasekhar has discussed that secular instability (which manifests 
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FIG. 3: A schematic diagram of the model problem. The curve given by Eq.JHJ passes through the origin for c = 0. For this 
critical trajectory O is a saddle point while C is a centre type fixed point. This figure differs from Fig. Q by a tilt of the axes. 



itself only if some dissipative mechanism is operative) arises in a Maclaurin spheroid, when the stresses derive from 
an ordinary viscosity which is defined in terms of a coefficient of kinematic viscosity (as the a parametrisation of 
Shakura and Sunyaev is for an accretion disc) , and when the effects arising from viscous dissipation are considered as 
small perturbations on the inviscid flow and taken into account in the first order |l9j . It is exactly in this spirit that 
we have prescribed the "quasi-viscous" approximation in the thin accretion disc. 

V. TEMPORAL EVOLUTION AS A SELECTION MECHANISM FOR TRANSONICITY 

In trying to make a time-dependent study, we have by now assured ourselves that a linearised perturbative analysis 
in real time cannot give us any conclusive insight about the accreting system showing any preference for any particular 
solution. If at all the stability analysis has indicated anything worth a serious consideration, it is that the presence of 
viscosity tends to destabilise the thin disc system. One may then say that any selection mechanism based on explicit 
time-dependence has to be non-perturbative and evolutionary in character. We now proceed to make a study of 
that, confining our discussion to only the inviscid disc, because we have seen that the steady solutions of the inviscid 
model are stable under small perturbations, and we may therefore expect stationary behaviour in the long-time limit. 
However, even for this simple inviscid system, with explicit time-dependence taken into consideration, the equations 
of a compressible flow cannot be integrated. Hence, to have any appreciation of the time-evolutionary selection of 
the critical solution, it should be instructive to consider an analogous model situation first. The model system being 
introduced here, describes the dynamics of the field y(x,t) as 




(33) 



whose static limit leads to 



dy y + 2x + x 2 
Ax y — x 



(34) 



and which, viewed as a dynamical system, is seen as 



-f- = y + 2x + x 2 



dr 



d:r 



(35) 




In the y — x space, the fixed points (x c , y c ) are to be found at (0, 0) and (—3, —3). A linear stability analysis of 
the fixed points in r space gives the eigenvalues A, by A 2 = 1 + 2(x c + 1). It is then easy to see that (0, 0) is a saddle 
point while (—3, —3) is a centre type point. 
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The integral curves are 



y 2 — 2xy — 2x 2 x 3 = c 



(36) 



and the trajectories passing through the saddle point are the ones with c = 0. The different possible trajectories 
(drawn without arrows associated with them) are shown in Fig. |3| The separatrices are shown as the heavy solid 
curve and the heavy dashed curve. The similarity between Fig. ^ and Fig. [21 is pretty much in evidence. 

To explore the temporal dynamics, and to obtain a solution to Ea. (|33[l . it would be necessary to apply the method 
of characteristics [2^ ■ This involves writing 



dt 



dx 



dy 



y — x y + 2x + x 2 



(37) 



The task is to find two constants c\ and from the above set and the general solution of Eq. would then be 
given by c\ = F{c2), where the function F is to be determined from the initial conditions. It is easy to see that one 
of the constants of integration is clearly the c of Ea. (|36[l . Hence, writing c\ = c and on using Ea. (|36|) in the first part 
of we get 



dt = ± 



dx 



^3x 2 + (2/3) x 3 + c 



(38) 



which solves the problem in principle. To put this in a usable form, the integration in Eq. l|38[l would have to be 
carried out. This cannot be done exactly. For small x (the most important region, since it is near the saddle), the 
x 3 term may be left out to a good approximation. Further, only the positive sign in the right hand side of Ea. H38|) 
is to be chosen by the physical argument that the system is to evolve through a positive range of t (time) values. 
Integration of Eq. (j38(l will then lead to the result 



(x + ^/x 2 + c/3) e" 



V3t 



C2 



(39) 



which will then make the solution of Eq. (|33|l look like 



y 2 - 2xy - 2x 2 



r 




(y - x)* 



— nr.3 



-1/3 1 



The evolution of the system is to be followed from the initial condition that y = for all x at t = 0. 
x 3 term again for small x, allows us to determine the form of the function F as F(z) — —3(2 — \fS)z 2 . 
consequently, becomes 



(40) 



Dropping the 
The solution, 



2xy - 2x 2 



r 



3 (2- 




(y-xY 



and as t — >oo, the steady solution that would be selected would be 



— x 3 
9 



-1\Tit 



(41) 



y 2 - 2xy - 2x 2 - -a; 3 = 



(42) 



which is actually the equation for the separatrices. It is worth stressing the remarkable feature of this result. The 
evolution started under conditions far removed from transonicity. In fact, it started as y = for all x (in the vicinity 
of the origin of coordinates) at t = 0. The evolution proceeded through a myriad of possible steady state solutions 
(all arguably stable under a linear stability analysis) and then in the stationary limit, selected the separatrices. This 
is a convincing demonstration that it is in principle possible for apparently non-realisable separatrices in the steady 
regime, to become eminently realisable physically, when the temporal evolution of the system is followed. 

To the extent that the model problem has been a good representative of the true physical situation, the whole 
argument for a time-dependent and non-perturbative method of selection developed above, may now be extended to 
the actual problem of thin disc accretion. For a steady accretion disc many previous works have upheld the case for 
transonicity Q, l2l| , although without explicitly addressing the issue of what special physical criterion may select the 
transonic solution to the exclusion of all other possible solutions. For the case of disc accretion on to black holes, Liang 
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and Thompson |2lJ make a clear point by saying that "the solution for the radial drift velocity of thin disk accretion 
onto black holes must be transonic, and is analogous to the critical solution in spherical Bondi accretion, except for 
the presence of angular momentum". For the inviscid disc at least, this is a most crucial analogy. In our chosen disc 
model, of course, the gravitational attractor driving the accretion process is not a black hole, but the question of 
the transonicity of the inflow solution is not expected to be too radically affected by this difference, especially in the 
vicinity of the far-off outer critical (saddle) point. However, what is important for us is the similarity with spherical 
symmetry that Liang and Thompson plj have indicated. The inviscid and thin disc model has rendered the related 
mathematical problem very conveniently one-dimensional, and has made the analysis very similar to a spherically 
symmetric flow. From this analogy it should be instructive for us at this point to derive some useful physical insight. 

Transonicity is a settled fact 0, Q in spherically symmetric accretion, and this has been so because at every spatial 
point in the flow, velocity evolves at a much greater rate through time, in comparison with density, and this therefore 
excludes the possibility of matter accumulating in regions close to the surface of the accretor. As a result gravity wins 
over pressure at small distances and the system is naturally driven towards selecting the transonic solution. This then 
raises the question of whether or not a likewise time-evolutionary mechanism should be at work for the selection of 
the critical inflow solution in an inviscid and thin accretion disc. We contend that it should be so. 

To have any appreciation of how the temporal evolution acts as a selection mechanism in disc accretion, it should 
be necessary to consider the dynamic equation for velocity evolution in the inviscid and thin accretion flow. This is 
given by 



dv dv , , 2 d P GM L 2 

8t +v m + w m + (43) 



It would be easy to see that to have a solution pass through the outer critical point (the saddle point), the temporal 
evolution should proceed in such a manner, that the inflow velocity would increase much faster in time than density 
at distances given by R ~ R C 2, where for sub-Keplerian flows, gravity going as R~ 2 , dominates the rotational effects, 
which go as R~ 3 . That the velocity starts decreasing in the region where R < R c x, is because of the resistance arising 
due to centrifugal effects in the fluid flow, which is quite unrelated to the density. The role of the pressure term 
(deriving from density) to resist the flow is confined to the subcritical solutions only. 

But what should be the key physical criterion guiding this dynamic selection of the transonic solution? We stress 
once again that the selection principle would very likely be the same as in the spherically symmetric case, where the 
transonic solution is chosen by dint of its corresponding to a configuration of the lowest possible energy |3j • To have 
an understanding of this, we consider Ea. i(43|) . with the density gradient being neglected as a working approximation. 
On large length scales, where the flow is highly subsonic, this is an especially effective approximation, and it allows 
us to treat the evolution of the velocity field to be largely independent of the density evolution. We then integrate 
the rest of the partial differential equation for velocity by the method of characteristics [2(J and get a solution that 
can be written down as 



v 2 GM L 2 



R 2R 2 



F 



vR Ct 

Ra + R(1 + v/C) \CR~ ~ R~q 



■ exp 



(44) 



in which C is an integration constant and Rq — GM/C 2 . The form of the function F is to be determined from the 
physically realistic initial condition v = at t = for all R, which will render F as 



GMS, | L 2 i 2 

+ 2(l-£Ro) a 



HO = -^r + - "L, 2 (45) 



If we examine the argument of F in Ea . ((44(1 and then study its long-time behaviour, we will see from Ea. 1)45(1 that for 
t — ► oo, the selected solution will correspond to 

v 2 GM L 2 , 
-2-jT + 2R 2 =° (46) 

We can now see that prior to the evolution, the system had no bulk motion and that the radial drift velocity was 
given flatly everywhere by v = 0. This, of course, gives the condition that initially the total specific mechanical 
energy of the system was zero. Then at t = both a gravitational mechanism is activated in this system and some 
angular momentum is imparted to it. This will induce a potential —GM/R everywhere, and at the same time start a 
rotational motion, respectively. The system will then start evolving in time, with the velocity v at each point in space 
evolving temporally according to Ea. 1)44(1. Finally the system will restore itself to a steady state in such a manner 
that the total specific mechanical energy at the end of the evolution (for t — > oo) will remain the same as at the 
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beginning (at t = 0), a condition that is given by Eo. l|46[) . whose left hand side gives the sum of the specific kinetic 
energy, the specific gravitational potential energy and the specific rotational energy. This sum is zero, and therefore, 
under the given initial condition, this must be the steady state corresponding to the minimum possible total specific 
energy of the system. Hence, this is the configuration that is dynamically and non-perturbatively selected. It is now 
conceivable that if the pressure term were to be taken into account, then the solution that would be dynamically 
selected, would be the one that would pass through the saddle point, since, as Bondi had analogously conjectured for 
spherically symmetric accretion |2| , this would be the one to satisfy the criterion of minimum energy. 

Interestingly enough, this question can also be addressed from a very different perspective. We have already 
derived an equation for a perturbation in the flow, as given by Ea. (|28|l . For the inviscid disc that we are studying, an 
intermediate step in arriving at this expression, with a = 0, is 



0_ 

dt 



fo V dt 



d_ 

dt 



df 

fo \dR 



d_ 

dR 



r„,2 



df 
fo V & 



d_ 

dR 



Jo 



dR 







all of whose terms can be ultimately rendered into a compact formulation that looks like |23 

d„ (l^duf) = 



(47) 



(48) 



in which we make the Greek indices run from to 1, with the identification that stands for t, and 1 stands for R. 
An inspection of the terms in the left hand side of Eq. I|47|l will then enable us to construct the symmetric matrix 



pi/ 



Vo 
fo 



1 



Vo 

vo v 2 -(3 2 c 2 Q _ 

Now the d'Alembertian for a scalar in curved space is given in terms of the metric g M „ by [25 

1 



(49) 



(50) 



with g^ u being the inverse of the matrix implied by g^ u . We use the equivalence that = A / = gg A " / , and therefore 
g = det (f Miy ), to immediately set down an effective metric for the propagation of an acoustic disturbance as 



=cff _ 



1 



Vo 



v v 2 -p 2 c 2 



(51) 



which can be shown to be entirely identical to the metric of a wave equation for a scalar field in curved space-time 22]. 
The inverse effective metric, g e Jl, can be easily obtained from Eq. l|51[l . and this will give v 2 = j3 2 c 2 as the horizon 
condition of an acoustic black hole [22|. Exploiting this close correspondence between the physics of supersonic 
acoustic flows and many features of black hole physics, we can argue that since infalling matter crosses the event 
horizon of a black hole maximally, i.e. at the greatest possible speed, by analogy the same thing may be said of 
the manner in which a rotating flow in an inviscid disc crosses an analogous acoustic horizon, given by the critical 
condition Vq — P 2 c 2 Q . This could only mean that the flow will cross the acoustic horizon transonically, since it is 
with the transonic flow that the greatest possible flow rate is associated. That we should be able to appreciate this 
fact through what is in essence a perturbative method, as given by Ea. H47|) . is quite remarkable, because conventional 
wisdom tells us that we should be quite unable to have any understanding of the special status of any inflow solution 
solely through a perturbative technique 

Having gained this insight into the preference for the transonic solution in an inviscid disc, it should also be quite 
obvious for us to realise that the introduction of any viscosity-dependent term in Eq. l|47|l will disrupt the precise 
symmetry of Ea. H48|) . With this disruption of the symmetry obtained from inviscid conditions, we should also lose 
all clear-cut vision of the way in which matter will cross the acoustic horizon. That this is what it should be in the 
presence of dissipation in the system, could also be understood from other considerations. The quasi- viscous disc being 
a dissipative system, i.e. energy being allowed to be drained away from this system, there should be no occasion for us 
to look for the selection of a particular solution, and a selection criterion thereof, on the basis of energy minimisation, 
as we have been able to conjecture for the idealised inviscid disc. This, of course, shall also affect the flow rate, and 
the whole system would be left without any well-defined criterion by which it could guide itself towards a steady 
state. Secondly, already having acquainted ourselves with the fact that the quasi-viscous disc is unstable on large 
length scales, we should expect no solution — transonic or otherwise — to be free of time-dependence. Therefore, a 
long-time evolution of the disc towards a stationary end is not something that we might hope for. 

Nevertheless, we are at least well aware of the fact that viscosity determines the distribution of matter in a viscous 
disc. This, of course, is intimately connected with the cumulative transfer of angular momentum to large length scales 
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of the disc. If we refer to Eg . (|15|) . we shall see that in the outer regions of the disc, the effective specific force is given 
to a first order in a by 



CM T 2 T 2 

+ (52) 



from which it is evident that on scales of R ~ Rl, the transport of angular momentum will give rise to an asymptotic 
constant non-zero force opposed to gravity. For the viscous disc this gives rise to a repulsive effect on large length 
scales vis-a-vis gravitational attraction. From this argument we go a step beyond and conjecture that a _1 / 3 i?L defines 
a limiting length scale for accretion that the outward transport of angular momentum imposes. This actually should 
be quite compatible with how viscosity redistributes an annulus of matter in a Keplerian flow around an accretor : 
the inner region drifting in because of dissipation, and consequently, through the conservation of angular momentum 
and its outward transport, making it necessary for the outer regions of the matter distribution to spread even further 
outwards 0, 0] . This state of affairs is qualitatively not altered in anyway for the quasi- viscous flow, except for the 
fact that with viscosity being very weak here, the outward transport of angular momentum can perceptibly cause an 
outward drift of matter only on very large scales. It is obvious that once a = 0, i.e. for the inviscid limit, this scale 
would be shifted to infinity. 



Acknowledgments 

This research has made use of NASA's Astrophysics Data System. One of the authors (AKR) gratefully acknowl- 
edges the support provided by the Council of Scientific and Industrial Research, Government of India, for a part of 
the time that was needed to carry out this work. The authors would also like to thank Dr. Tapas Kumar Das, Prof. 
Rajaram Nityananda and Prof. Paul J. Wiita for some helpful comments. 



[1] S. K. Chakrabarti, Theory of Transonic Astrophysical Flows (World Scientific, Singapore, 1990). 
[2] H. Bondi, Monthly Notices of the Royal Astronomical Society 112, 195 (1952). 
[3] A. R. Garlick, Astronomy and Astrophysics 73, 171 (1979). 

[4] M. A. Abramowicz and W. H. Zurek, The Astrophysical Journal 246, 314 (1981). 
[5] R. Yang and M. Kafatos, Astronomy and Astrophysics 295, 238 (1995). 

[6] D. Molteni, H. Sponholz, and S. K. Chakrabarti, The Astrophysical Journal 457, 805 (1996). 

[7] P. J. Wiita, Accretion Disks around Black Holes in Black Holes, Gravitational Radiation, and the Universe (Kluwer, 
Dordrecht, 1998). 

[8] T. K. Das, J. K. Pendharkar, and S. Mitra, The Astrophysical Journal 592, 1078 (2003). 

[9] J. Frank, A. King, and D. Raine, Accretion Power in Astrophysics (Cambridge University Press, Cambridge, 1992). 
[10] N. I. Shakura and R. A. Sunyaev, Astronomy and Astrophysics 24, 337 (1973). 
[11] R. Matsumoto, S. Kato, J. Fukue, and A. T. Okazaki, Publ. Astron. Soc. Japan 36, 71 (1984). 
[12] M. Afshordi and B. Paczyhski, The Astrophysical Journal 592, 354 (2003). 

[13] D. W. Jordan and P. Smith, Nonlinear Ordinary Differential Equations (Clarendon Press, Oxford, 1977). 

[14] A. K. Ray, Monthly Notices of the Royal Astronomical Society 344, 83 (2003). 

[15] R. Narayan and I. Yi, The Astrophysical Journal 428, L13 (1994). 

[16] J. E. Pringle, Annual Review of Astronomy and Astrophysics 19, 137 (1981). 

[17] S. Kato, F. Honma, and R. Matsumoto, Monthly Notices of the Royal Astronomical Society 231, 37 (1988). 
[18] X. Chen and R. Taam, The Astrophysical Journal 412, 254 (1993). 

[19] S. Chandrasekhar, Ellipsoidal Figures of Equilibrium (Dover Publications, New York, 1987). 

[20] L. Debnath, Nonlinear Partial Differential Equations for Scientists and Engineers (Birkhauser, Boston, 1997). 

[21] E. P. T. Liang and K. A. Thompson, The Astrophysical Journal 240, 271 (1980). 

[22] M. Visser, Classical and Quantum Gravity 15, 1767 (1998). 



